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Abstract 

We consider an infinite network of globally-coupled phase oscillators in which the natural fre- 
quencies of the oscillators are drawn from a symmetric bimodal distribution. We demonstrate that 
macroscopic chaos can occur in this system when the coupling strength varies periodically in time. 
We identify period-doubling cascades to chaos, attractor crises, and horseshoe dynamics for the 
macroscopic mean field. Based on recent work that clarified the bifurcation structure of the static 
bimodal Kuramoto system, we qualitatively describe the mechanism for the generation of such 
complicated behavior in the time varying case. 
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In nature and in many practical applications, it is not uncommon to observe the emer- 
gence of coherent macroscopic behavior in large populations of interacting rhythmic units 
despite noise and the presence of heterogeneity in the population. For systems of globally- 
coupled phase oscillators, collective synchrony and simple macroscopic oscillations have been 
extensively studied [l-6]. Large networks of more complicated (e.g., chaotic) oscillators can 
also exhibit these behaviors [7], but can also display a collective chaotic state js]. It is not 
clear, however, if heterogeneous networks constructed of simple phase oscillators that are 
not independently capable of producing chaos can exhibit complex macroscopic behavior 
(such as chaos) in the thermodynamic limit of large system size. In the current work, we use 
a recently-developed mean-field analysis method and demonstrate that chaos can exist in 
the macroscopic mean field for a heterogeneous network of globally coupled phase oscillators 
with a bimodal frequency distribution and time-periodic coupling. We propose a qualitative 
mechanism for how this arises, and we identify period-doubling scenarios, attractor crises, 
and Smale horseshoe dynamics. 



I. INTRODUCTION 



A remarkably successful and analytically tractable model for describing the spontaneous 
onset of coherence in large populations of phase oscillators was introduced by Kuramoto 

n n n n 

in 1975 |2|, |3[, and many subsequent extensions have been formulated p, |6[. Typically, 
Kuramoto-like models share the following three fundamental characteristics: (i) The indi- 
vidual rhythmic units within the network are simple phase oscillators. When isolated, the 
temporal evolution of the oscillators is given by 9i = Ui, % = 1, ■ ■ ■ ,N, where Ui is the intrin- 
sic natural frequency of the i-th oscillator, and the number of oscillators in the network (N) 
is assumed to be large, (ii) The collection of natural frequencies is assumed to be distributed 
according to a time- invariant function g(u). (iii) The coupling among oscillators is assumed 
to be all-to-all (i.e. global), so that each oscillator influences, and is influenced by, all others 
in the network. 

The original Kuramoto model assumes a smooth unimodal distribution function g(oj) 
which is symmetric about a mean frequency uq and which monotonically approaches zero as 
\u— u>o\ — > 00. When uncoupled, the network is incoherent, and each oscillator drifts in phase 
with respect to the others. But, when the coupling strength is sufficiently strong (and N is 
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sufficiently large), a coherent state emerges through a continuous phase transition at a critical 
coupling strength K*. A macroscopic domain of phase- locked oscillators begins to form and 
coherence grows as the coupling strength K continues to increase. Traditionally, this phase 
transition is quantified by a complex order parameter z (to be defined in Sec. (Ill AD ) which 
describes the macroscopic mean field. This has magnitude zero in the incoherent state and 
becomes nonzero for K > K*. 

Previous efforts analyzing the emergence of coherence have mainly focused on the loss of 
stability of the incoherent state. A recent breakthrough introduced by Ott and Antonsen 
(9), 1(3] has allowed us to move beyond this approach, and there have been many develop- 



ments in understanding the Kuramoto system and its extensions [HHl5|. The OA method 
identifies a low- dimensional invariant manifold to which the dynamics of the macroscopic 
mean field of a large (N — > 00) heterogeneous network of globally-coupled phase oscillators 
is attracted. Remarkably, low- dimensional equations of motion for the mean-field behavior 
on this manifold can be derived. For example, the governing equation for the order param- 
eter of the original Kuramoto problem (on the manifold) can be written as a single complex 
nonlinear ordinary differential equation. Recent analyses of a number of Kuramoto-type 
models using this method have revealed a rich set of possible nonlinear dynamical states 



16| : limit cycles, chimeras, quasi-periodic states, multistability, and standing waves, as well 



as various local and global bifurcation scenarios for the mean field, including saddle node, 
transcritical, Hopf, saddle- node infinite-period (SNIPER), and homoclinic bifurcations. 

But can a network of simple phase oscillators result in more complicated dynamical be- 
havior, such as chaos, in the infinite-iV limit? Kuramoto himself speculated that "within 
the framework of the phase model. . . no chaotic dynamics at the collective level seems pos- 



sible" 



17j |. For globally-coupled complex Ginzburg-Landau-type oscillators, clustering and 



macroscopic chaos have been observed 17H20|. However, it has been argued that this re- 



quires a degree of freedom in the oscillator amplitude [13, ll8| , and thus the model diverges 
from the simple phase oscillator description. Chaotic behavior has also been reported in 
coupled map lattices js|, but the individual elements of such systems are typically nonlinear 
maps that are separately capable of producing chaos when uncoupled. There is numerical 
evidence suggesting that chaos exists in finite populations of a number of Kuramoto-type 



phase models 



21 



22|, but this appears to disappear as iV — > 00. Chaos in the order 



rameter has also been reported for a resistively-loaded Josephson junction array 



23 



pa- 



24]. 



3 



However, the analytic model in 24j corresponds to the non-generic case of a Kuramoto-like 
system with a homogeneous (delta function) distribution of natural frequencies. This is 
notable because it has been shown that the dynamics for the classical Kuramoto system 
9] and the resistively loaded Josephson junction network jiol . \\\ cannot posses a chaotic 
attractor in the infinite- iV limit in the more generic situation with a spread in the oscillator 
natural frequencies. 

The question is then: can attracting macroscopic chaos or other complicated dynamics 
exist in a heterogeneous network of simple phase oscillators in the thermodynamic limit 
(N oo)? If so, what might the minimum requirements be? Here, we provide a partial 
answer by extending the bimodal Kuramoto model, which we analyzed previously 11], 
to the case in which the global coupling parameter is a periodic function of time. We 
derive equations of motion for order parameters which are valid for N — > oo and long 
times, and we use these to firmly establish the existence of chaos in the macroscopic mean 
field by identifying period-doubling cascades to chaos, attractor crises, and Smale horseshoe 
dynamics. 



II. FORMULATION 



n n 

The bimodal Kuramoto system has been investigated by Kuramoto [3|], Crawford [4|, and 
251 ] . For the case of a bimodal natural frequency distribution g{oS) consisting 



others 



m 



of the sum of two offset but otherwise identical Cauchy-Lorentz distributions, a complete 



bifurcation diagram was recently reported 



ll|. (The same authors also found that when 



the sum of two offset Gaussians was used for g(u), the resulting bifurcation diagram was 
qualitatively the same.) 

Most related work assumes that the connectivity of the network of interest is static. 
However, in natural networks (i.e., physical, biological, social, etc.), communication among 
the individual components is often time- dependent (see 26( and the references therein). The 
coupling strength, the type of coupling, and the connection topology may not necessarily 
remain constant. Here, we generalize the bimodal Kuramoto system to include the situation 
in which the global coupling parameter varies periodically in time. 
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A. Our Model 



We consider the following system: 

K(t) N 

0* = + £ sin^- - ft) (1) 

3=1 

for % = 1, 2, . . . , N, in which the coupling strength is assumed to be a periodic function given 
by 

K(t) = K + Asm (^ytj . (2) 

The strength and the period of the periodic variation are given by A and r respectively 
These two parameters will serve as the bifurcation parameters in the current study. 

The natural frequency Ui for each oscillator is randomly drawn from a normalized bimodal 
distribution function given by the sum of two Cauchy-Lorentz distributions, 

9{U) = 27 ((w- J)2 + A* + (c + coWa 2 ) ' (3) 
where A characterizes the half-width of the individual peaks and ±uq are their center fre- 
quencies (see Fig. [TJ). Without loss of generality, we choose g(u) to be symmetric about 
u = 0, as one can always pick an appropriate rotating frame in which this holds. However, 
for g(u) to be effectively bimodal, one must choose co > A/y^ so that the central dip of 
the bimodal distribution is convex. A more general study of the Kuramoto system with 
an asymmetric bimodal distribution will be reported elsewhere j^, and a case in which a 
bimodal form of q(oj ) that cannot be written as the sum of two even unimodal distributions 
was analyzed in [15]. 



B. OA reduction 



In the infinite- N limit, we can describe the phase oscillators within our network at a 
particular time t by a continuous distribution function F(9,u,t), such that F(9,u,t)d9du 
gives the fraction of phase oscillators with phases between 9 and 9 + d9 and with natural 
frequencies between u and u+du. Since the number of oscillators is assumed to be conserved 
at all times, we have 

/OO /*27T 
/ F(9,u,t)d9du=l, 
-oo JO 
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and since g{ui) is assumed to be constant in time, the marginal distribution density 

F(6,u,t) d9 = g(u) 

re evolution of Kuramoto-type systems in terms 



10(). In particular, the time evolution of F must 



is as well. Many authors have analyzed t 
of this dilution function (e,.,jfi 3 
satisfy the continuity equation, 

where the phase velocity vg is given by the continuum version of Eq. [I], 

2tt 



(4) 



vo(0,u,t)=u + K(t) / F(9,co,t)sm{9' -9)d9. 



(5) 

The macroscopic mean field is described by a complex order parameter originally introduced 
by Kuramoto [3j, 



z(t) = pe 



iip 



oo p2n 



F(9,u,t)e ie dddu. 



(6) 



'-oo JO 

Geometrically, the order parameter describes the centroid of all the phasors e td within the 
network. When the network is incoherent, z(t) is zero; it becomes nonzero when coherence 
emerges. Our goal is to describe the potentially complicated dynamical states for this 
macroscopic variable when we vary the system parameters A and r. 
The phase velocity vg can be re- written in terms of z(t) as follows: 



v e {0, u, t) = G + i [H(t)e- W - H*(t)e l6 ] 



(7) 



where H(t) = K{t)z{t) and G — to. Note that due to the global coupling, H(t) does not 
explicitly depend on the individual phase 9 except through the mean field z(t). In fact, one 
can take Eq. (j7|) as the starting point for a generalization to a larger class of Kuramoto-type 
systems for which the OA reduction method is applicable. The method applies as long as 
the network vector field can be put into the single harmonic form given by Eq. where 
neither H nor G depend on 9 explicitly. 

Eqs. (jSJ), and <^7§ form the governing equations for the distribution function F(9, u, t). 



Following Ref. 



9], we expand F(9,u,t) in a Fourier series according to the ansatz 



F(9,u,t) 



2tt 



l + ^a n (w,t)e 

n=l 



i nO 



+ ^a* n (w,t)e" 

n=l 



ind 



(8) 
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where a(u,t) is a yet-to-be-determined function independent of 9 and which has modulus 
less than one (so as to guarantee convergence of the Fourier series). It is important to note 
that this ansatz defines a submanifold within the infinite-dimensional space of all possible 
distribution functions F. By directly substituting the ansatz, Eq. OH]), into the continuity 
equation, Eq. (j3J, one obtains 



+ iua + i [Ha 2 - H*]\ ^n« n - 1 e m9 + c.c. = 0, 

^ - ' n=l 



(9) 



where c.c. denotes the complex conjugate of the expression on the left- hand- side of the 
equation. Then, since Yl^Li na n ~ l e m9 = jrz^my. 7^ 0, the factor within the parentheses 
must vanish identically if the ansatz given by Eq. flSJ) is to be a valid solution to the continuity 
equation. This then gives 

3cv 1 

— + icoa + - \Ha 2 -H*} =0, (10) 
at 2 L J 

which holds for each value of u. Note that H will in general depend on the macroscopic 
mean field z(t), which can be written 

/OO /»27T ^OO 
/ F(6,u,t)e w d6du = / a*(u,t)g(w)du. (11) 
-oo J J —oo 

Therefore, the time evolution of a(u,t), and thus the distribution function F(8,u,t), can 
be obtained by integrating the integro-differential equation given by Eqs. f lTUj) . (ITT]) and 
the mean field function H(z) [2j|. Note that although the use of the ansatz has effectively 
collapsed the time evolution of the infinite number of Fourier modes in 6 into one evolution 
equation for a, Eqs. (ITOl and ( ITT]) remain an infinite-dimensional system since we have a 
continuous set of natural frequencies u. As we will show in our specific examples below, 
Eqs. ( ITO]) and ( ITTj) reduce further to only a small number of ordinary differential equations. 



C. No Chaos in the Classic Kuramoto Model 

The original Kuramoto model can be put into the form of Eq. (JJJ) if we define H(t) = 
Kz(t), G = u, and let u be chosen randomly from a unimodal Cauchy-Lorentz distribution. 
We then have 

v e (6,u,t) = u + ^[ze- ie -z*e ie ] . 
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Then, following the formalism described above, the amplitude a(w,t) evolves in time ac- 
cording to 

da K r 2 -, 

— + iuja + — \ za - z = 0, 

at 2 L J 

with the mean field z(t) being given by Eq. (ITT]) . Following Ref. |9|, we further assume that 
a(u,t) can be analytically continued into the lower complex w-plane and that the initial 
conditions satisfy (i) |a(u;,0)| < 1, and (ii) |a(w,0)| — > as $s(u)) —> oo. Then, using a 
semi-circular contour in the lower complex w-plane to perform the integral in Eq. ( TTTj) . and 
taking the radius of the contour to infinity, we can express z(t) in terms of the residue of 
the contour integral, i.e., z(t) = a*(uo — iA,t). Therefore, at large time, the macroscopic 
mean field for the classical Kuramoto system evolves according to a single complex ordinary 
differential equation: 

^ = -(A + ^ + f (z-z*z 2 ). (12) 

Since this is a two-dimensional system, no complicated behavior such as chaos is possible 
for the classical Kuramoto system with a heterogeneous distribution of natural frequencies 
in the thermodynamic limit. This result does not contradict the numerical results obtained 



in Refs. 



21 



22], since those examples considered finite-size networks. In fact, numerical 

n 

evidence in Ref. [21J indicated that the largest Lyapunov exponent in the finite-size network 
tends toward zero as the system size increases. 



D. Bimodal Kuramoto Model with Time-varying Coupling 

The same functions H(t) and G used above apply to the bimodal case with time- 
dependent coupling. The difference is the form of the distribution function g(u). Following 
the same procedure outlined above leads to 

^ + iuja + ^- \za 2 -z*} =0, (13) 
at 2 L J 

where the coupling K(t) is now time-dependent according to Eq. (T5]) and the mean field z(t) 
is given by Eq. (JTTj) . Since explicit time dependence appears only in H(t) = K(t)z(t) = 
(K + Asm(27rt/T)z(t), the mathematical analysis is similar to the time-independent case 
reported in [llj. We briefly summarize the procedure and the relevant results. 

The bimodal Lorentzian distribution has two simple poles in the lower w-complex plane, 
u = ±cuo — iA, and direct integration of Eq. (TiTi) gives z(t) = 21 M+ Z2 M with 21,2 (t) = 

8 



a*(±u>o — iA,t). These two sub-order parameters represent the two populations of phase 



osci 
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lators whose natural frequencies cluster around u = +uq and uj = — uiq, respectively 
291 ] . Substituting this expression for z(t) into the differential equation for a, we arrive 
at the following pair of nonautonomous equations for Z\(t) and z 2 (t): 

(14) 



<f = -(A + iu ) Zl + ^1{ Z1 + Z2 - (zl + z *) z f)} 



^ = -(A - iu )z 2 + m [zi + Z2 - ( Z * _ z * )z i)] . 

We expect the asymptotic behavior of these two complex sub-order parameters to be sym- 
metric except for a relative phase difference ip such that ^ = e 1 ^ and \z 2 \ — — P (for 
details, see jlO] and [llj]). Substituting this symmetry condition into Eq. (JTi]) . we obtain 
equations that describe the long-time macroscopic mean-field behavior for the time- varying 
bimodal Kuramoto system: 



p = -Ap + ^p(l-p 2 )(l + cos^) 
j) = 2u - ^(l + p 2 )sinV. 



(15) 



Note that the OA reduction of the time-static network in [llj has the same mathematical 
form as this, except that here, K(t) is a function of time. The full bifurcation structure 
of these equations with K(t) = constant was analyzed in llj, and the main features are 
shown in Fig. [2j For convenience, we refer to this system as the "static" system, and we 
write "time- varying system" to describe the case when K(t) varies in time according to 
Eq. ©. 

An important result of Ref. [llj, which analyzed the static version of Eqs. (Tl5|) . is that 
two particular combinations of parameters, Auq/K and 4A/K, reveal the full dynamical 
structure of the system. Figure [2] shows the parameter space described by these expressions. 
The incoherent state is stable for parameter values at the top of this figure. Its stability 
can be lost either via a transcritical bifurcation (TC, semicircular boundary) or by a Hopf 
bifurcation (HB, half-line). A saddle-node (SN) bifurcation occurs in the vicinity of where 
the TC and the HB line meet (but away from the incoherent state in state space), and 
thus there is an approximately triangular region of multistability in which an attracting 
partially-synchronized state coexists with the attracting incoherent state. The multistable 
region extends below the Hopf curve to the homoclinic (HC) bifurcation curve until the SN 
curve joins it and becomes a SNIPER curve. 
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Eqs. ( JT5j) describe the macroscopic mean field of the full network (Eqs. (JTK2J) ) with iV — > oo 
once the dynamics has converged onto the manifold defined by the ansatz (Eq. [8]). Ref. 
showed that convergence onto this ansatz manifold occurs as long as there is heterogenity 
in the natural frequencies of the oscillators. We note that Eqs. ffTBT) describe the dynamics 
on the ansatz manifold for both the static and the time-varying system, since the the time- 
dependent coupling that we introduce in Eq. ([2]) carries through in the reduction procedure. 
In particular, K(t) does not produce excursions transverse to the ansatz manifold. However, 
under circumstances to be examined below, it produces complex trajectories within the 
ansatz manifold. Thus, we expect these complex trajectories to be observable in the full 
network. 

III. MACROSCOPIC DYNAMICS FROM A TIME- VARYING BIMODAL NET- 
WORK 

Below, we use Eqs. ( !T5|) to analyze the asymptotic dynamics of the macroscopic mean 
field (p, ip) with the time-dependent coupling strength given by Eq. fl2]). To understand these 
dynamics, it is useful to consider the relationship between the attractors of the static network 
and those of the time- varying system. Recall that in the latter, the coupling strength varies 
periodically according to K(t) = K + Asm (2nt/r). Thus, as K(t) varies (with fixed values 
of Wo and A), the parameters in Figure [2] vary and sweep through one of the regions indicated 
by the diagonal green lines. The attractors of the static system become "moving targets" in 
the time-varying system. We will refer to these "moving targets" as pseudo-static attractors. 
The resulting behavior of the time- varying system depends on the amplitude A and frequency 
1/t of the time- varying coupling. We consider the five cases labeled (a) through (e) marked 
by the red points in Figure [2J which indicate the static parameter values (uq, A) for the 
average value of K(t) in each case. We set Kq = 4 throughout. See the figure caption for 
the remaining parameter values. 

A. The Incoherent State and Periodic and Quasiperiodic Solutions 

In cases (a)-(c), we set r = 5 and A = 0.3. Note that the amplitude A is sufficiently weak 
such that the time-dependent system parameters remain fully within each respective region 
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of stability relative to the static network. 

For case (a), both the static system and the time-dependent system possess a fixed attract- 
ing equilibrium at p = corresponding to the incoherent state. Although the eigenvalues 
of this equilibrium oscillate in the time-varying system, it remains an attracting state, and 
hence the long-time behavior is not affected by the time variation of the coupling. 

At the parameter values for case (b), the static system is attracted to an equilibrium 
that corresponds to a partially-synchronized state with p > 0. In contrast to case (a), this 
equilibrium is destroyed by the time variation of the coupling. In its place, the time- varying 
system exhibits a periodic orbit with period r. This arises as the time-varying system's 
trajectory follows the "moving target", the pseudo-static equilibrium. The link between 
equilibria of the static network and periodic orbits of the time-varying network (away from 
bifurcations) can be made rigorous by applying standard averaging methods (see, e.g., [30 1 



or 



ail). 



Upon crossing the Hopf bifurcation line from above, the incoherent state of the static 
network loses stability, and an attracting periodic orbit emerges. This periodic orbit is the 
only attractor in case (c) for the static netowrk. For the time- varying system, the frequency 
of variation 1/r is typically not commensurate with the (effective) intrinsic frequency of 
oscillation of this periodic orbit. Thus, the asymptotic state for the time- varying network 
in this case is generally a quasiperiodic orbit on a torus |32|. We also observed frequency 
locking behavior in this case as the amplitude A was varied (not shown) 33]. 

We numerically confirmed that the time- dependent reduced system, Eq. ( fl5|) . and the 
full system, Eq. both exhibit these local dynamical features. In the latter case we 
used an ensemble of 500, 000 oscillators. As expected, the weakly time- varying networks 
appropriately "follow" the expected behavior of the static Kq networks. Fig. [3] compares 
the asymptotic behavior of the macroscopic mean field for the reduced and the full systems 
in cases (a)-(c). In case (a), both systems converge to the incoherent state, but the full 
system exhibits fluctuations on the order of 1/yN due to finite-system-size effects. The 
same effects can be seen in case (b), in which the full-system limit cycle is "thickened" by a 
similar amount. In case (c), the apparent thickness of the ring-like attractor is much larger 
and reflects deterministic quasiperiodic dynamics as predicted by the reduced analysis (see 
inset). 

These results are not surprising in light of the averaging theorem 34| and the results 
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of Ref. 



10] . But, we emphasize that the dynamics induced by the periodic variation in 



K(f) remain within the ansatz manifold, and we have confirmed that the behavior of the 
time- dependent reduced equations is evident in the full system. Of greater interest is the 
behavior of the time-varying system when the parameters sweep across bifurcations. We 
examine this in the next section. 



B. Period Doubling Cascades, Chaos, and a Crisis 

Case (d) is different from the previous cases in that the parameters sweep across multiple 
bifurcations as well as a small region of multistability. We note that as K = Kq+A sin(27rt/r) 
varies in the range [Kq — A,Kq + A] with Kq = 4 and A = 1, the parameters effectively 
sweep back and forth along the green line across this region for case (d) in Fig. [2j To clarify 
the situation, we show in Figure H] (A) a one-dimensional bifurcation diagram, calculated 
using the static system with K = Kq ± A, in which the magnitude of the asymptotic 
macroscopic mean-field p is plotted versus K. The left panel shows the static asymptotic 
structures at each fixed parameter value K. The incoherent state is present throughout, and 
is attracting for K < 3.2. At K — 3.2, the Hopf bifurcation is encountered, and a periodic 
orbit (whose minimum and maximum p values are plotted) emerges. At K = 3.935, a saddle- 
node bifurcation creates two equilibria, and the unstable one collides with and destroys the 
periodic orbit in a homoclinic bifurcation at K = 3.953. For larger values of K, only the 
equilibria exist. Fig. H](B) shows all three static structures in state space for K = 3.94, when 
they co-exist (X = pcosip and Y = psinip are plotted). 

The dynamics of the time-varying system can be understood in terms of transitory at- 
traction to the pseudo-static attractors as K(t) varies. For slow sweeping (r small), there is 
ample time for the state of the time- varying system to approach the pseudo-static attractors 
as they drift in time. Thus the system either stays close to the pseudo-static attracting 
equilibrium or approximates the pseudo-static periodic orbit, depending on which is present 
at any given time. For extremely slow sweeping, this results in behavior similar to a bursting 
neuron, in which periods of quiescence alternate with "bursts" of near-periodic excursions 
in phase space (see the movie in Fig. (fTTj) in the Appendix). For very fast sweeping, the 
pseudo-static attractors come and go too quickly to affect the dynamics, and the time- varying 
system's behavior approaches that of the time-averaged system, i.e., the static system with 
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K = (K(t)) = Kq (see the movies in Figs. (11 2p and (|T3|) . Appendix) [26l. |35|. For intermedi- 
ate sweeping frequencies, however, transitory attraction to moving pseudo-static attractors 
can dominate the dynamics of the time- varying system and complicated dynamics can arise, 
as we show below. 

Using the reduced system, Eq. (1151) . we first investigate the time- varying network behav- 
ior by fixing r = 5 (corresponding to an intermediate sweeping frequency) and gradually 
increasing the amplitude A of the periodic coupling strength variation from zero. Accord- 
ingly, we sweep across the pseudo-static attractors shown in Figure IHA), from — A to A, for 
increasing values of A. We investigate changes in behavior associated with r later in this 
section. 

For A = 0, the time-varying system reduces to the static system with (u , A, K ) = 
(1.29,0.8,4). In this case, the static system is attracted to an equilibrium with p / 0. For 
small values of A > 0, the macroscopic mean-field exhibits a small limit cycle (libration) 
circling near the static equilibrium (not shown, but similar to Figs. [3]B). As A increases, 
a saddle-node bifurcation occurs at A = 0.441, creating a new stable limit cycle with full 
rotation in phase space. This limit cycle, shown in Fig. EJ^A) for A = 0.55, is 1:1 phase- 
locked to the drive period r. As A increases further, a period-doubling occurs at A = 0.558 
(Fig.[5]^B)), and then again at A = 0.600 (Fig.[5](C)). This period-doubling cascade continues, 
and by following the sequence through period sixteen and using the Feigenbaum constant 
36L 1371 ] . we estimate the infinite-period accumulation point to be at A^ = 0.615. Beyond 
this, chaos is found. Fig. E](D) shows the chaotic attractor obtained for A = 0.65 (see the 
movie in Fig. ( JT41) . Appendix). 

It is easier to visualize this bifurcation sequence using a Poincare surface of section with 
ip — ijjQ chosen appropriately, so that limit cycles appear as fixed points. We chose trajectory 
crossings through ip = 1.1 radians with ij) > and obtained the bifurcation diagram shown 
in Figure El The cascade described above is clearly visible in the range 0.41 < A < 0.62. 



Antimonotonicity is also evident 38]: chaos is destroyed through a sequence of reverse 
period-doubling bifurcations in the range 0.82 > A > 0.982, returning to a full-rotation, 1:1 
limit cycle. We have also observed several other limit cycles with different locking ratios 
created through similar saddle-node bifurcations. These orbits follow their own period- 
doubling cascades into chaos (e.g., near A = 0.5075 (1:4) and A = 0.529 (1:3); note that 
these are small and are not clearly apparent in the figure). Note also that the small limit 
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cycle (libration) created in the weak time-variation regime, visible on the left near p = 0.68, 
terminates at a saddle node bifurcation at A = 0.545. 

The lower panel of Figure [5] shows the two largest Lyapunov exponents, calculated using 
the reduced system (Eq. [T5"|) . Regions in which one is positive are clearly evident, as are 
reg ions in which neither is positive, e.g., in one of the infinitely- many periodic windows 



39|,|40|. 

Focusing on the main 1:1 phase- locked bifurcation branch, we see that the chaotic bands 
merge in the usual way and become a one-piece chaotic attractor near A = 0.627. Then, 
at A c = 0.641, the chaotic attractor explodes in size through an interior crisis. This occurs 
when the attractor touches the stable manifold of a coexisting unstable period-three orbit. 
Figure [7(A) shows a stroboscopic map in (p, tp) (points sampled every r time units) showing 
the smaller chaotic attractor when A < A c . For A > A c , a typical orbit remains near the 
core smaller attractor most of the time, but intermittently, bursts occur during which the 
trajectory visits a more expanded region of state space (see Fig. [7J (B)). A sequence of 
iterates (numbered red dots) during one of the these bursts is shown in panel (a) as the 
orbit transiently visits the mediating unstable period-three orbit. Interestingly, this crisis 
results in a loss of synchronization between the relative phase difference ip and the phase of 
the coupling strength variation. Before the crisis (Fig. [7J(A)), the attractor is restricted to a 
limited range of ip in the stroboscopic map (note the axes), indicating that the macroscopic 
dynamics is phase-locked to the external periodic modulation. After the crisis (Fig. [7J(B)), 
the attractor's phase ip ranges throughout the entire to 2ir range. Thus, the macroscopic 
mean- field slips in phase with respect to the coupling modulation 4ll |. 

We explicitly demonstrate the existence of a chaotic set in the macroscopic mean field 
using the stroboscopic map described above. Setting A = 0.65, we choose the set of initial 
conditions (p(nr), ip(nr)) indicated by the curved trapezoidal region in Fig. [HI The next 
iterate (p((n + l)r),ip((n + l)r)) of this set under the stroboscopic map (equivalent to 
integrating these initial conditions forward by r time units) forms the elongated and curved 
region shown in the figure. These intersect in the manner of a Smale horseshoe, thus 
indicating that the reduced equations for the time-varying network contain a chaotic set 
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All of the above examples had the period r of coupling variation set to 5. Continuing 
with the parameters for case (d), we fix A = 0.65 and examine the effect of changing r. 
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The results appear in Fig. [9j which shows a bifurcation diagram of p versus r obtained via 
Poincare surface of section as above (again with ipo = 1.1 radians subject to ip > 0). For 
intermediate values of r, the macroscopic mean-field once again exhibits period-doubling 
cascades, antimonoticity, crises, and periodic windows. At the extremes, the expected be- 
havior discussed earlier is evident. For slow coupling variation (6.67 < t ~ 16.0), the system 
exhibits a simple periodic orbit, in which the trajectory follows the pseudo-static attracting 
equilibrium and limit cycle well; this appears as the fixed point on the right of Fig. [9j For 
very fast coupling variation (r < 4.1), the fixed point on the left of Fig. [9] corresponds to a 
very small periodic orbit (libration) that approximates the equilibrium of the static system 
with K = (K(t)) = Kq. This can be seen in the movies of Figs. [T21 and [T51 in the Appendix, 
in which the static system's equilibrium is marked with an "X" . 

Finally, we briefly consider case (e), in which we return to varying A with r = 5. Here, 
the parameter sweep only crosses the SNIPER bifurcation (see Fig. [2]). Figure [TPl shows that 
period-doubling cascades to chaos, crises, and regions with positive Lyapunov exponents are 
once again evident, although in a smaller range of A. 

IV. DISCUSSION 

We have demonstrated that macroscopic chaos occurs in the thermodynamic limit of the 
bimodal Kuramoto system with periodic time-variation of the coupling strength. This is 
interesting because it is known that an infinite but autonomous network of phase oscillators 
with a symmetric bimodal natural frequency distribution cannot exhibit chaos. This follows 
from analysis based on the OA reduction procedure, which shows that the latter system is 
effectively two-dimensional at the level of the mean field dynamics. Thus, chaos is impossible. 
However, by introducing a periodic time- variation of the coupling strength, a third dimension 
is added to the system, and under proper circumstances, chaos is found. 

To clarify the circumstances under which macroscopic chaos is possible in the time- varying 
system, we have developed the concept of pseudo-static attractors that act as "moving 
targets" which influence the trajectory of the macroscopic mean field. These pseudo-static 
attractors are the attractors of the static system for any given fixed coupling strength K. As 
K varies in time, these become pseudo-attractors that move about in state space depending 
on the frequency (or period) of variation. As a result, the trajectory of the time-varying 
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system may or may not be able to "keep up" with the motion of the pseudo-attractors. We 
find that interesting complicated dynamics arise when the frequency of variation is such that 
the trajectory of the macroscopic mean field is essentially frustrated in that it is not able 
to ever reach the moving pseudo-attractors. In this situation, the trajectory is dominated 
by what we call transitory excursions. The frequencies for which complicated dynamics 
arise should, in some sense, be commensurate with the rates of attraction to the pseudo- 
static attractors. In the examples that we have described here, we observed period-doubling 
cascades to chaos, crises, horseshoes, and other complex nonlinear dynamical features. 

Our understanding of this mechanism is based on detailed knowledge of the asymptotic 
structures and bifurcations of the macroscopic mean field present in the static bimodal Ku- 
ramoto system, published earlier [ll|. Our observations suggest the following conjecture. In 
order to generate macroscopic chaos on the antsatz manifold, the static system must exhibit 
the following dynamical features: (1) Topological Choice: Within the range of paramet- 
ric variation (in our case, Kq ± A), the static system must possess at least two attracting 
macroscopic structures that are separated in state space. The separation is necessary in 
order to permit substantial transitory excursions. In our examples (case (d) and (e)), the 
static system has an attracting equilibrium and/or an attracting limit cycle. (2) Switch- 
ing: Within the range of parameter variation, there must exist bifurcations of the static 
system that create and destroy these attracting macroscopic structures. In our examples, 
the attracting equilibrium is created/destroyed through a saddle-node bifurcation, and the 
attracting limit cycle is created/destroyed through a homoclinic bifurcation. Furthermore, 
the parametric variation should include a range in which these macroscopic attractors exist 
independently. This is necessary in order to cause transitory excursions. For case (d) (see 
Fig. H] (A)), the equilibrium is the only attractor in the static system for K > 3.953, and 
the limit cycle is the only attractor in the static system for K G [3.2, 3.935]. Although both 
of these attractors coexist in a very small interval for case (d), there is no such coexistence 
in case (e), demonstrating that multistability is not necessary for complex behavior to arise 
in the time- varying system (Fig. [TOl) . 

Finally, we note that this general mechanism need not be restricted to non-autonomous 
systems. We suspect that a similar mechanism could arise in any system that has a dynamic 
interplay between fast and slow dynamics under the conditions described above. Indeed, we 
have found complicated behavior, including chaos, in an autonomous model of a neuron 
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featuring fast spiking dynamics which are modulated by slow autonomous ion concentration 
dynamics 43|. We leave this for future investigation. 



V. APPENDIX 

The following movies, available online, help to visualize the dynamics described in the 
main text. In all movies, case (d) is considered (ojq = 0.8, A = 1.29) with A = 0.65. 
The pseudo-static attractors are shown in black (compare Fig. and the trajectory of the 
time- varying system is shown in cyan, being traced by the red dot. 

Figure [11] shows an extreme case of slow variation (r = 50), and is reminiscent of a 
bursting neuron. Figures [121 and [T3l show the fast (r = 1) and very fast (r = 0.1) variation 
case, illustrating the time- varying system's approach to the equilibrium of the static system 
with K = (K(t)) = K . The location of the latter is marked with a black "X". Finally, 
the chaotic trajectory obtained for r = 5 is shown in Fig. [TU 
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FIG. 1: Bimodal distribution of natural frequencies, g(w), as a sum of two Lorentzians. 
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FIG. 2: (Color online) Phase diagram for the static bimodal Kuramoto network with coupling K = 
Kq. The incoherent state is stable in the upper region, and black curves denote bifurcations that 
lead to coherent collective states (TC=transcritical, HB=Hopf, SN=saddle-node, HC=homoclinic, 
SNIPER=saddle-node-infinite-period). In the time-varying system, Kit) = Kq + A sin(27rf/r), and 
the parameters sweep along the short diagonal green lines. The lettered red points indicate cases of 
interest for which (wo,A) are as follows: (a) (2.0,1.5); (b) (0.7,1.35); (c) (2.5,0.8); (d) (1.29,0.8); 
(e) (1.4, 0.65). We set Kq = 4 and we have A = 0.3 for cases (a)-(c) and A = 1 for case (d). 
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FIG. 3: Behavior of the reduced (left) and full (right) time-varying equations for cases (a)-(c) with 
Kq = 4, r = 5, and A = 0.3. A: Case (a), the persistence of the incoherent state. The fluctuations 
in the full system are due to finite-system-size effects. B: Case (b), a stable fixed point in the static 
system becomes a limit cycle (libration) in the weakly time- varying system. C: case (c), a limit 
cycle in the static system becomes a quasi-periodic state in the weakly time-varying system. 
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FIG. 4: (Color online) A: Bifurcation diagram showing the static asymptotic structures past which 
we sweep. Solid lines are stable equilibria; dashed lines are unstable equilibria; lines with symbols 
represent the maxima (circles) and minima (squares) of limit cycles. B: State space diagram 
showing the limit cycle and the equilibria (circles; solid for stable and open for unstable) for 
K=3.94, when these coexist. X = pcostp, Y = psintp, and in both panels, A = 0.8 and loq = 1.29. 
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FIG. 5: Period doubling and chaos in the time- varying network for case (d). The panels show Y = 
psinip versus X = p cos V for A=0.55 (A), 0.56 (B), 0.61 (C), and 0.65 (D), with u = 1.29, A = 0.8, 
and Kq = 4. 
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FIG. 6: (Color online) A: Bifurcation diagram, obtained using a Poincare surface of section, showing 
the magnitude p of the macroscopic mean field versus the amplitude A of coupling variation, in 
case (d). B: Plot of the two largest Lyapunov exponents versus A. Both panels were computed 
using the reduced system, Eq. (fT5j) . 
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FIG. 7: (Color online) a) The chaotic attractor of the macroscopic mean- field slightly above and 
below the interior crisis value at A c = 0.641. For these diagrams, a stroboscopic surface of section, 
(p((n+ 1)t), ip((n + 1)t)) vs. (p(nr), ip(riT)), was used. A: the attractor for A < A c . Superimposed 
on this is a trajectory segment for A > A c showing escape via the mediating period-three orbit. 
B: the expanded attractor for A > A c . The pre-crisis attractor is overlaid in blue. 
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FIG. 9: The bifurcation diagram for the magnitude p of the macroscopic mean-field versus the 
period r of the coupling time variation. A Poincare surface of section is used. A = 0.65, ojq = 0.8, 
and A = 1.29. 
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FIG. 10: (Color online) Bifurcation diagram and Lyapunov exponents for case (e); compare Figj6j 
u = 1.4, A = 0.65. 
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FIG. 11: Very slow parameter variation with r = 50.0 (enhanced online). 
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FIG. 12: Fast parameter variation with r = 1.0 (enhanced online). 
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FIG. 13: Very fast parameter variation with r 


= 0.1 (enhanced online). 
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FIG. 14: Chaos: intermediate parameter variation with r = 5.0 (enhanced online). 
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